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1.  INTRODUCTION 


The  objective  of  this  contract  is  to  develop,  incorporate,  test,  and  validate  new 
algorithms  for  Nascap-2k  that  are  needed  to  self-consistently  compute  plasma  transport 
and  to  model  electromagnetic  radiation  in  the  near  to  mid  field  from  VLF  (3  kHz  to 
30  kHz)  antennas.  The  plasma  flow  models  can  be  used  to  address  various  plasma 
engineering  concerns  including  surface  discharges  due  to  meteoroid  impact  and 
spacecraft  contamination  due  to  electric  propulsion  plasma  plume  effects.  The  goal  of 
this  effort  is  to  provide  a  plasma  engineering  capability  to  the  spacecraft  community. 

During  the  third  and  fourth  year  of  this  contract,  new  database  and  memory  manager 
software  was  developed  for  Nascap-2k,  support  was  provided  for  the  DSX  program,  and 
an  algorithm  was  developed  for  the  computation  of  surface  currents  on  a  satellite  that  is 
acting  as  an  antenna. 

I.l.  Nascap-2k 

During  this  period,  the  primary  improvement  made  to  Nascap-2k  was  the 
development  of  a  new  database  and  memory  manager  system.  This  is  discussed  in 
Section  1.2. 

In  March  2008,  we  provided  an  interim  delivery  of  Nascap-2k  (version  3.1.3)  to 
AFRL.  In  January  2009,  before  incorporation  of  the  new  database  and  memory  manager 
software,  we  released  Nascap-2k  3.2  to  AFRL  and  the  NASA  SEE  Program,  which 
distributes  the  code.  Prior  to  release,  we  tested  the  code  on  our  standard  suite  of  problems 
on  Windows  and  Linux  in  order  to  confirm  that  recent  changes  had  not  introduced  any 
problems. 

In  order  to  test  and  develop  the  multiprocessor  capabilities  of  Nascap-2k,  we  ported 
the  code  to  a  4  processor  Apple  Macintosh  Pro  with  the  MacOS  X.5  operating  system 
with  Intel  C++  and  Fortran  compilers.  Nascap-2k  is  now  fully  compatible  with  the 
Windows,  Red  hat  LINUX,  and  MacOS  X  environments. 

We  wrote  N2kScriptRunner,  a  C++  code  that  runs  a  Nascap-2k  script  outside  of  the 
Java  user  interface.  Using  N2kScriptRunner,  we  compiled,  linked,  and  ran  Nascap-2k 
using  the  OpenMP  compiler  commands  for  multiprocessor  operations. 

We  implemented  an  improvement  to  the  algorithm  used  to  display  nodal  extensive 
quantities  such  as  the  charge.  The  approach  used  is  described  in  Section  2. 

In  addition  to  the  above,  we  made  a  number  of  small  changes  to  Nascap-2k.  The  most 
notable  are  the  following: 

•  We  modified  the  Results3D  tab  of  the  Nascap-2k  user  interface  to  allow  a  user  to 
request  that  a  parameter  be  displayed  on  multiple  planes  at  once. 


•  We  added  to  the  Results  tab  the  ability  to  view  the  surface  number  of  the  surface 
with  the  minimum  and  maximum  value  of  the  displayed  quantity. 

•  We  added  to  the  Results  3D  tab  the  ability  to  view  the  components  of  the  current. 

•  We  added  the  Freja  environments  as  default  auroral  environments  on  the 
Environment  tab. 

•  We  made  the  specification  of  an  insulating  surface  as  “fixed  potential”  work 
properly. 

•  We  removed  the  use  of  schema  files. 

•  We  added  a  check  to  make  sure  that  the  particle  file  filename  is  no  more  than  20 
characters  (the  most  that  can  be  read  by  the  keyword  input  routines  used  by 
Tracker). 

•  We  added  a  check  to  let  the  user  know  that  there  is  a  missing  conductor  number. 

1.2.  New  Database  and  Memory  Manager 

During  the  fourth  year  of  this  contract,  we  focused  on  designing,  building,  and 
implementing  into  Nascap-2k  the  new  database  and  memory  management  software, 
N2kDB.  N2kDB  is  complete  and  has  been  incorporated  into  Nascap-2k  as  an  initial 
implementation.  In  order  to  simplify  the  transition,  the  initial  implementation  uses 
wrappers  between  the  original  database  commands  and  the  new  ones.  Thus,  in  the  short 
term,  some  of  the  constraints  of  the  old  database  remain.  Over  the  next  few  months,  all  of 
the  old  database  commands  will  be  replaced  with  new  ones. 

We  began  by  writing  a  Software  Design  Document.  This  document  describes  the 
software  design,  test  procedures,  and  software  standards.  The  document  was  revised 
during  development  as  the  design  matured.  This  document  has  subsequently  been 
reorganized  into  N2kDB  Software  Documentation.  This  document  appears  as  Volume  II 
of  this  report. 

We  developed  N2kDB  Test,  a  database  testbed  code.  N2kDB  Test  has  the  same 
structure  as  Nascap-2k  and  uses  the  database  in  the  same  manner.  It  is  a  small 
manageable  code  that  we  used  for  testing  the  database  software  during  development  and 
for  developing  the  proper  implementation  of  the  new  commands.  We  anticipate  using  it 
intermittently  over  the  next  few  months;  however,  we  do  not  plan  to  incorporate  it  into 
the  released  code. 

We  developed  two  versions  of  N2kDBTool  (console-based  and  with  a  Java 
interface),  a  stand  alone  program  that  reads  and  writes  Nascap-2k  database  files.  This 
program  has  proved  invaluable  during  development  of  N2kDB  and  promises  to  be  useful 
in  future  Nascap-2k  development.  N2kDBTool  will  be  included  with  future  releases  of 
Nascap-2k. 

N2kDB  also  includes  two  testing  programs  (msiotest  and  dmtest)  that  verify  the 
behavior  of  the  database  code  itself  They  verify  that  the  requested  operation  is  performed 
correctly.  (N2kDB  Test  verifies  that  the  appropriate  operation  was  requested.)  These 
codes  will  be  maintained  with  N2kDB,  but  will  not  be  distributed  to  users. 
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N2kDB  satisfies  all  the  requirements  specified  in  the  Software  requirements 
document.  A  review  of  how  these  requirements  are  satisfied  is  given  in  an  appendix  to 
N2kDB  Software  Documentation. 

Before  incorporation  of  N2kDB  into  Nascap-2k,  we  released  Nascap-2k  3.2. 

1.3.  DSX  Calculations 

We  used  Nascap-2k  to  perform  a  set  of  self-consistent  calculations  of  the  plasma 
response  to  a  high  voltage  square  wave  VLF  antenna  using  the  capabilities  implemented 
earlier  and  described  in  the  previous  interim  report  for  this  contract.  These  calculations 
are  included  in  a  paper  prepared  for  the  Spacecraft  Charging  Technology  Conference  and 
are  the  first  full  test  of  all  of  the  new  capabilities.  This  paper  (M.J.  Mandell,  V.A.  Davis, 
D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Nascap-2k  Self-consistent  Simulations  of  a  VLF 
Plasma  Antenna,  Spacecraft  Charging  Technology  Conference,  Biarritz,  France,  June 
2007)  is  included  in  Appendix  A.  We  learned  that  the  inclusion  of  a  representation  of  the 
thermal  distribution  at  the  boundary  can  influence  the  current  collected  by  the  antenna 
arms. 

1.4.  Surface  Currents 

We  are  developing  a  technique  for  the  evaluation  of  surface  currents  for  DSX  and 
prototyped  it  in  Java.  A  description  of  the  algorithm  and  the  prototype  implementation  is 
in  Section  3. 

1.5.  DSX  Program 

We  supported  the  DSX  program  by  assuming  an  active  role  in  the  Y-boom  high 
voltage  isolation  design.  To  that  end.  Dr.  Myron  J.  Mandell  participated  in  the  following 
teleconferences  and  meetings. 

•  Weekly  teleconferences  with  Y  Antenna  supplier  (L’Garde)  from  November  6, 
2007  to  through  Februray  13,  2008. 

•  A  visit  to  the  L’Garde  facility  to  discuss  the  design  issues  on  November  7,  2007. 

•  Y  Antenna  Isolation  teleconference  on  December  19,  2007. 

•  Weekly  teleconferences  with  the  alternate  supplier  of  Y  Antenna  (ATK)  from 
January  15  to  February  25,  2008. 

•  ATK  Isolation  face-to-face  at  Houston  Airport  (lAH)  on  February  13,  2008. 

•  DSX  Alternative  Y  Antenna  CDR,  Goleta,  California  on  February  26  by 
teleconference. 

•  DSX  Y  Antenna  CDR  at  L’Garde,  Tustin,  California  on  February  27. 
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Dr.  Myron  Mandell  attended  and  made  a  presentation  at  the  Workshop  on  The 
Remediation  of  Enhanced  Radiation  Belts  in  Lake  Arrowhead,  California  on  March  3-6, 
2008. 

Dr.  Mandell  also  attended  the  DSX  System  CDR,  Breckenridge,  Colorado,  May  6-8, 
2008. 

Dr.  Mandell  attended  the  MURI  Review  and  RBR  Workshop  at  Stanford  in  Palo 
Alto,  California,  Februrary  18-19,  2009. 

1.6.  Contract 

The  scientists  and  other  researchers  who  contributed  to  this  work  are  as  follows:  Dr. 
Myron  J.  Mandell,  Dr.  Victoria  A.  Davis,  Ms.  Barbara  M.  Gardner,  Ms.  Katherine 
Wilcox,  Ms.  Alisa  J.  Ward,  and  Mr.  Nicholas  R.  Baker. 

This  contract  is  a  follow-on  to  work  performed  under  earlier  contracts:  F 19628-9 1-C- 
0187,  Space  System-Environment  Interactions  Investigation;  F19628-93-C-0050, 
Modeling  and  Post  Mission  Data  Analysis;  F19628-89-C-0032,  Analysis  of  Dynamical 
Plasma  Interactions  with  High  Voltage  Spacecraft;  and  F19628-98-C-0074,  Spacecraft 
Potential  Control.  NASA  supported  related  work  under  contracts  NAS8-98220  and 
NAS8-02028. 

1.7.  Publications 

The  following  publications  were  supported  in  total  or  in  part  by  this  contract  during 
the  third  and  fourth  years. 

M.J.  Mandell,  V.A.  Davis,  D.L.  Cooke,  A.T.  Wheelock,  C.J.  Roth,  Nascap-2k  Self- 
consistent  Simulations  of  a  VLF  Plasma  Antenna,  Spacecraft  Charging  Technology 
Conference,  Biarritz,  France,  June  2007. 

M.J.  Mandell,  V.A.  Davis,  E.J.  Pencil,  M.J.  Patterson,  H.K.  McEwen,  J.E.  Foster,  J.S 
Snyder,  Modeling  the  NEXT  Multi -Thruster  Array  Test  with  Nascap-2k,  Spacecraft 
Charging  Technology  Conference,  Biarritz,  France,  June  2007.  (Original  research 
supported  by  NASA.) 

M.J.  Mandell,  V.A.  Davis,  E.J.  Pencil,  M.J.  Patterson,  H.K.  McEwen,  J.E.  Foster,  J.S 
Snyder,  Modeling  the  NEXT  Multi-Thruster  Array  Test  with  Nascap-2k,  IEEE 
Transactions  on  Plasma  Science,  36,  p.  2309,  2008.  (Original  research  supported  by 
NASA.) 

M.J.  Mandell,  V.A.  Davis,  A.  Wheelock,  D.L.  Cooke,  Modeling  Space  Weather  Effects 
using  Nascap-2k,  GOMACTech-08,  Las  Vegas,  NV,  March  2008. 
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2.  NODAL  CHARGE  DENSITY 


Using  an  earlier  version  of  Nascap-2k,  attempts  were  made  to  plot  the  value  of  the 
Nodal  Charge  Density  along  a  line.  The  plots  had  non-physical  structures,  particularly 
near  grid  boundaries.  Also,  the  Nodal  Charge  Density  on  a  plane  as  displayed  on  the  3-D 
Results  tab  showed  non-physical  values  along  grid  boundaries.  The  reason  for  the  non¬ 
physical  values  is  that  the  Nodal  Charge  Density  as  displayed  is  actually  the  charge  on 
each  node  as  stored  in  the  database  and  used  in  the  calculations  (an  extensive  quantity) 
converted  to  the  charge  density  (an  intensive  quantity)  for  the  purpose  of  display. 
Extensive  quantities  exist  only  on  nodes,  and  therefore  it  does  not  make  sense  to  discuss 
their  value  at  an  arbitrary  point.  The  correct  way  to  convert  an  extensive  quantity  to  an 
intensive  quantity  is  essentially  to  multiply  by  the  inverse  of  the  volume  matrix. 
Unfortunately,  while  straightforward,  implementing  this  is  quite  complex.  With  a 
moderate  level  of  effort  we  were  able  to  implement  an  improvement  to  the  display  of 
nodal  extensive  quantities.  The  approach  used  is  as  follows: 

•  In  order  to  usefully  display  an  extensive  quantity,  it  is  necessary  to  convert  it  to  an 
intensive  quantity. 

•  Previously,  the  conversion  was  done  by  dividing  the  value  at  each  node  by  the  cube 
of  the  mesh  unit.  This  would  be  correct,  if  there  were  only  one  grid  (except  at  the 
outer  boundary  and  near  the  object). 

•  We  added  two  grid  interface  operations  to  the  cut  plane  display  routines.  One  grid 
interface  operation  is  before  and  one  after  the  division  of  each  node  by  mesh  volume. 
In  the  first  grid  interface  operation,  the  extensive  quantity  on  nodes  on  grid 
boundaries  is  transferred  from  the  iimer  grid  to  the  outer  grid.  In  the  second,  the 
intensive  quantity  (resulting  from  the  division  by  volume)  is  shared  from  the  outer 
grid  to  the  iimer  grid.  No  adjustment  is  made  to  account  for  the  fact  that  the 
appropriate  volume  for  grid  interface  nodes  and  those  near  the  object  is  not  the  mesh 
unit  cubed. 

The  nodal  charge  density,  derived  from  the  nodal  charge  using  the  old  and  new 
approaches  for  display  is  shown  on  the  following  figures  for  the  moving  sphere  problem 
provided  as  an  example.  Most  of  the  artificial  structures  on  the  grid  boundaries  have  been 
eliminated. 
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Figure  1.  Nodal  Charge  Density  on  X  =  0  plane  using  original  and  revised  algorithm 
for  display  of  extensive  quantity. 


Figure  2.  Nodal  Charge  Density  on  Z  =  0  plane  using  original  and  revised  algorithm 
for  display  of  extensive  quantity. 


6 


3.  SURFACE  CURRENT  ALGORITHM  AND  PROTOTYPE 
IMPLEMENTATION 

3.1.  Objective 

The  objective  of  the  project  is  to  calculate  surface  currents  on  a  satellite  that  is  acting 
as  an  antenna.  The  satellite  has  two  or  more  antenna  elements  (represented  by  conductor 
numbers)  that  are  biased  in  a  programmed  way  by  an  amplifier,  and  typically  will  have 
passive  elements  as  well.  The  amplifier  injects  current  at  a  specified  location  on  each 
antenna  element  in  order  to  maintain  the  programmed  voltage.  The  injected  currents, 
along  with  plasma  currents  resulting  fi'om  the  voltage  application,  flow  on  the  spacecraft 
surface.  These  currents  are  the  source  of  radiated  electromagnetic  fields. 

The  plan  is  to  add  to  Nascap-2k  the  capability  to  compute  these  radiated 
electromagnetic  fields.  Nascap-2k  divides  spacecraft  surfaces  into  surface  elements.  The 
fields  could  be  computed  by  adding  the  contributions  fi'om  the  current  elements  (vectors) 
for  each  surface  element.  The  current  element  for  each  surface  would  be  the  average 
surface  current  density  times  the  element  area  and  located  at  the  element  centroid. 

We  have  developed  an  algorithm  to  calculate  these  surface  currents  and  are  creating  a 
prototype  implementation  in  Java  prior  to  incorporation  into  Nascap-2k. 

3.2.  Algorithm 

Each  antenna  or  passive  element  is  represented  by  a  number  of  surface  elements  with 
a  common  conductor  number.  The  amplifier  provides  the  charge  needed  (calculated 
based  on  the  change  in  electric  field  and  the  plasma  currents)  to  achieve  the  specified 
potential  change  of  the  conductor  during  the  time  step.  This  charge  is  injected  into  a 
specified  surface  element  and  can  then  flow  across  the  element  edges  into  neighboring 
surface  elements.  For  surface  elements  other  than  the  injection  element,  the  change  in  the 
charge  on  the  element  must  equal  the  sum  of  the  plasma  currents  and  the  surface  currents 
crossing  its  edges,  as  shown  in  Figure  3.  The  current  crossing  each  edge  is  iteratively 
computed  (as  described  below).  The  element  surface  current  (the  average  current  on  the 
element)  is  the  current  that  best  matches  the  edge  currents.  This  process  is  performed 
separately  for  each  conductor.  Current  is  not  allowed  to  flow  across  an  edge  separating 
two  different  conductors. 


7 


Totalchangc  in  charge  Plasma  Currents  Surface  Currents 


Figure  3.  The  change  in  charge  on  a  surface  element  is  made  up  of  plasma  currents 
and  surface  currents. 

Figure  4  shows  the  steps  in  this  process.  The  injection  current  for  each  conductor  is 
determined  by  summing  over  all  the  surface  elements  of  the  conductor.  The  edge  currents 
are  initialized  by  each  element  contributing  the  appropriate  portion  of  its  required  rate  of 
change  in  charge.  The  edge  currents  are  then  adjusted  until  the  net  current  to  each  surface 
element  is  as  required.  At  the  begirming  of  the  “Augment  Edge  Currents” step,  we  know 
the  discrepancy,  Aleiem,  between  the  amount  of  current  leaving  the  element  across  its 
edges  and  the  amount  of  current  required.  The  current  crossing  each  edge  of  the  element 
is  augmented  by 


AI 


edge 


=  rAI 


elem 


L 


where  i  is  the  length  of  the  edge,  L  is  the  total  length  of  edges  bounding  the  element,  and 
r  is  a  relaxation  factor  currently  set  to  0.8.  Each  non-bounding  edge  also  receives  a 
contribution  from  its  neighboring  element.  In  our  experience  thus  far,  the  process 
converges  robustly. 
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Figure  4.  Iterative  calculation  of  edge  currents  to  obtain  surface  currents. 
This  algorithm  requires  the  following  quantities: 

1 .  The  surface  elements,  including  specification  of  the  injection  element  for  each 
conductor. 

2.  The  rate  of  change  of  external  electric  field  for  each  element. 

3.  The  plasma  current  to  each  surface  element. 


4.  The  change  in  charge  to  be  supplied  by  surface  current  is 


gpdE 

dt 


less  the  plasma 


current. 


3.3.  Prototype  Implementation 


The  prototype  is  pure  Java,  most  of  which  is  derived  from  existing  Nascap-2k 
interface  coding.  Most  of  the  work  is  done  by  the  SurfaceCurrents  class,  with  inputs 
supplied  by  the  two  auxiliary  classes  DEbyDT  and  PlasmaCurrents.  The  code  relies 
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heavily  on  existing  ObjectToolkit  classes,  notably  Mesh,  Element,  and  their 
dependencies.  Non-intrusive  modifications  have  been  made  to  Element  (e.g.,  to  store  the 
surface  currents  and  to  register  the  injection  surface  element)  and  to  Mesh  (e.g.,  for 
graphic  display  of  surface  currents). 

The  code  reads  in  the  object  mesh  and  separates  it  into  separate  meshes  (by 
conductor).  At  completion  the  surface  current  density  is  entered  into  each  element  of  the 
original  mesh  (as  a  VectorS).  Current  elements,  which  could  be  summed  to  calculate  the 
vector  potential,  are  given  by  the  surface  current  density  times  the  element  area. 

AL  =  rAI ,  -  (L 

edge  clein 

where  €  is  the  total  length  of  edges  bounding  the  element,  r  is  the  rate  of  change  of 
external  electric  field  for  each  element  is  computed  by  a  class  DEbyDT.  Presently 
DEbyDT  reads  the  .POTS  file  from  a  previous  Nascap-2k  calculation  using  the  recently 
developed  N2kDBTool.dll,  and  differences  the  stored  electric  field  at  adjacent  timesteps. 
The  use  of  actual  time  step  times  is  not  yet  implemented. 

The  plasma  current  to  each  surface  element  is  computed  by  the  class  PlasmaCurrents. 
Presently,  the  plasma  current  is  calculated  by  a  clsEnvironment  using  the  surface  element 
potential. 

3.4.  Graphics 

The  surface  current  (if  it  exists)  for  each  element  is  displayed  using  a  JavaSD 
primitive  Cone  object.  The  height  of  the  Cone  is  proportional  to  the  square  root  of  the 
surface  current  density  and  to  the  square  root  of  the  surface  element  area.  As  output,  the 
code  converts  the  entire  object  to  a  primitive  and  writes  it  out  along  with  the  element 
surface  currents.  The  primitive  can  then  be  displayed  by  Object  Toolkit.  The  same 
display  code  should  work  in  Nascap-2k  if  the  surface  elements  appear  in  the  mesh  that  is 
plotted. 

3.5.  Example 

For  development,  we  used  a  test  object  consisting  of  two  adjacent  cubes  (conductors 
1  and  2)  as  shown  in  Figure  5.  The  center  surface  elements  on  either  end  (shown  as  Gold) 
are  flagged  as  injection  surface  elements  (to  which  the  bias  currents  are  supplied).  During 

dE 

initial  development,  hard-wired  —  values  were  applied,  and  alternate  injection  locations 

dt 

were  used  to  demonstrate  that  the  process  converges  to  sensible  results  for  asymmetric 
cases. 

For  this  example,  a  Nascap-2k  run  was  performed  starting  with  conductor  1  at  OV  and 
conductor  2  at  -lOOOV,  and  the  bias  varied  sinusoidally  at  10  kHz  for  20  timesteps  of  5 
microseconds  each  (making  one  full  period).  The  surface  currents  application  was  then 
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run  for  two  cases,  and  the  output  saved  at  each  timestep.  We  then  recorded  the  time 
development  of  surface  current  at  surface  elements  78  and  31,  as  shown  in  Figure  5. 


The  two  cases  were:  (a)  a  small  amount  of  plasma  current  1 80°  out  of  phase  with  the 

E  dE 

applied  voltage;  and  (b)  plasma  current  comparable  to  the  —  current.  For  case  (a)  we 

dt 


expect  the  surface  currents  for  the  two  surface  elements  to  be  nearly  identical,  as  current 
is  injected  into  one  conductor  and  leaves  the  other  conductor  in  a  sinusoidal  pattern,  90° 
out  of  phase  with  the  applied  voltage.  For  case  (b)  we  expect  the  two  surface  currents  to 
be  different  and  phase  shifted,  as  the  plasma  current  is  always  outwardly  directed  (from 
the  center  of  the  object  toward  the  injection  points),  so  that  it  adds  to  one  while 
subtracting  from  the  other. 


Figure  6  shows  the  surface  currents  for  the  two  selected  surface  elements  in  the  low 
plasma  density  case(a).  The  small  amount  of  plasma  current  causes  the  two  curves  to  be 
slightly  separated,  and  results  in  a  phase  shift  of  about  15  degrees,  as  indicated  by  the 
accompanying  sine  wave. 


Figure  5.  Object  used  for  surface  current  test  and  development. 
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Figure  6.  Surface  current  for  selected  surface  elements  in  the  low  plasma  density 
case. 

Figure  7  shows  the  surface  currents  for  the  two  selected  surface  elements  in  the  high 
plasma  density  case(b).  The  significant  amount  of  plasma  current  causes  the  two  curves 
to  be  substantially  separated,  and  results  in  a  phase  shift  of  about  45  degrees.  The  glitches 
seen  near  50  microseconds  for  surface  element  78  and  near  100  microseconds  for  surface 
element  3 1  results  from  electron  plasma  current  collected  when  the  surface  elements  are 
near  zero  potential. 
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Figure  7.  Surface  current  for  selected  surface  elements  in  the  high  plasma  density 
case. 
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Figure  8  shows  surface  currents  for  the  high  density  case  at  timestep  7,  near  the 
dE 

maximum  in  — .  At  this  time  all  surface  currents  are  in  the  same  direction,  and  would  be 
dt 

of  similar  magnitude  for  the  two  conductors  if  the  plasma  current  were  not  significant. 

dE 

Figure  9  shows  the  surface  currents  for  the  same  case  near  a  null  in  — .  The  surface 
^  dt 


currents  are  now  due  solely  to  the  plasma  currents,  and  flow  in  opposite  directions  on  the 
two  conductors. 


Figure  8.  Surface  currents  for  the  high  plasma  density  case  at  timestep  7. 
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Figure  9.  Surface  currents  for  the  high  plasma  density  case  at  timestep  12. 
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APPENDIX  A:  yVA5C4P-2/irSELF-CONSISTENT  SIMULATIONS  OF  A  VLF 
PLASMA  ANTENNA 

The  following  paper  was  presented  at  the  Spacecraft  Charging  Technology 
Conference,  Biarritz,  France  in  June  2007. 
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ABSTRACT:  We  simulate  the  plasma  response  to  high  voltage  square 
wave  VLF  antenna  in  Medium  Earth  Orbit  plasma  with  Nascap-2k.  The 
plasma  is  modeled  with  a  hybrid  Particle-in-cell  (PIC)  approach  with  PIC 
ions  and  fluid  barometric  electron  densities.  The  plasma  response, 
collected  ion  currents,  and  chassis  floating  potential  are  computed  self- 
consistently  with  a  near-square-wave  bias  applied  to  the  antennas. 
Particle  injection  and  splitting  are  used  to  replenish  the  plasma  depleted  at 
the  boundary,  represent  the  thermal  distribution,  and  maintain 
appropriately  sized  macroparticles.  Therefore,  current  limitation  due  to  the 
thermal  distribution  of  ions  and  the  resulting  angular  momentum  barrier 
are  included.  Above  the  ion  plasma  frequency  the  plasma  current  lags  the 
voltage  by  about  10^,  while  below  the  ion  plasma  frequency  the  current 
leads  the  voltage  by  about  7° 


1  -  INTRODUCTION 

The  response  of  a  plasma  to  very  low  frequency  (VLF)  (3  kHz  to  20  kHz)  high-voltage  antennas  at 
orbital  altitudes  of  1000  to  10,000  kilometers  has  been  a  subject  of  scientific  interest  for  many 
decades.  ’  ’  Plasma  waves  from  VLF  antennas  with  such  frequencies  (“whistler”  waves)  are 
thought  to  interact  with  MeV  radiation  belt  electrons.  As  this  antenna  frequency  is  less  than  both  the 
electron  plasma  frequency  and  the  electron  gyrofrequency  (both  nearly  300  kHz  for  a  plasma 
density  of  10^  m‘^  and  a  magnetic  field  of  0.1  gauss),  only  certain  modes  can  propagate  as 
electromagnetic  waves,  and  the  near  field  is  dominated  by  electrostatic  effects.  Although  a 
comprehensive  self-consistent  electromagnetic-electrostatic  simulation  would  be  the  desired  goal, 
there  are  many  computational  challenges  to  be  overcome,  so  we  begin  with  electrostatic  simulations 
in  order  to  sort  out  the  dominant  plasmadynamic  effects. 
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A  VLF  dipole  antenna  has  two  elements,  each  several  centimeters  in  diameter  and  many  meters 
long.  Due  to  the  ease  of  electron  collection  by  positive  objects,  the  positive  element  is  at  near  zero 
or  small  positive  potential  relative  to  the  ambient  ionosphere,  while  nearly  the  full  applied  potential 
appears  on  the  negative  element.  Because  the  frequencies  of  interest  are  comparable  to  the  ion 
plasma  frequency,  the  sheath  structure  is  at  some  intermediate  state  between  the  “ion  matrix”  or 
“frozen  ion”  limit  (which  assumes  the  ions  are  stationary  and  contribute  ambient  ion  density  to  the 
space  charge)  and  the  equilibrium  space  charge  limit  (in  which  the  ions  assume  a  steady-state  space 
charge  limited  distribution  of  charge  and  current).  Thus,  calculation  of  the  sheath  structure  and  of 
the  ion  collection  by  the  antenna  requires  dynamic  (specifically,  particle-in-cell,  PIC)  treatment,  at 
least  for  the  ions.  Nascap-2k  can  be  used  to  perform  all  four  simulations  of  interest:  (1)  equilibrium 
space  charge  sheath;  (2)  “frozen  ion”  sheath;  (3)  dynamic  PIC  ions  with  fluid  (Boltzmann  or 
barometric)  electrons  (Hybrid  PIC);  (4)  dynamic  PIC  ions  and  electrons  (Full  PIC).  Previously  we 
reported  on  one  and  three  dimensional  Full  PIC  calculations'^  and  Hybrid  PIC  calculations  with 
prescribed  potentials  for  a  higher  density  plasma.^ 

This  paper  presents  antenna  simulations  performed  with  Nascap-2k.  The  plasma  is  modeled  using 
the  hybrid  PIC  approach  with  PIC  ions  and  fluid  barometric  electron  densities.  The  plasma 
response,  collected  ion  currents,  and  chassis  floating  potential  are  computed  self-consistently  with  a 
near-square-wave  bias  applied  to  the  antennas.  Current  limitation  due  to  the  thermal  distribution  of 
ions  and  the  resulting  angular  momentum  barrier  are  included.  The  different  plasma  response  above 
and  below  the  ion  plasma  frequency  are  shown. 

In  order  to  perform  these  calculations,  Nascap-2k''s  PIC  capability  has  been  expanded  to  include 
boundary  injection  and  particle  splitting.  The  boundary  injection  replaces  particles  that  are  either 
collected  or  leave  the  problem  in  long-running  calculations.  The  particle  splitting  allows  for 
modeling  of  the  effects  of  the  thermal  distribution  of  velocities,  as  well  as  reducing  particle  size 
when  passing  into  regions  of  finer  resolution.  The  algorithms  used  for  these  capabilities  along  with 
their  accuracy  and  limitations  are  discussed. 


2  -  PARTICLE  IN  CELL  TECHNIQUES 

Nascap-2k  has  primarily  been  used  to  model  quasi-static  phenomena.  However,  there  are  a  large 
number  of  physical  processes  of  interest  whose  timescales  require  a  dynamic  approach  such  as  a 
particle-in-cell  (PIC)  technique.  Examples  of  such  processes  are  breakdown  phenomena,  plasma 
kinetics,  and  sheath  structure  about  surfaces  with  potentials  that  change  on  a  timescale  comparable 
to  the  time  it  takes  an  ion  to  cross  the  sheath.  PIC  techniques  can  also  be  used  to  address  problems 
in  which  anal)^ic  representations  of  the  environmental  currents  are  inadequate,  such  as  in  a 
spacecraft  wake  or  in  a  cavity.  A  steady-state  PIC  technique,  in  which  the  ion  space  charge  density 
is  computed  from  macroparticles  tracked  from  the  boundary  of  the  computational  space  until  they 
are  collected  or  exit  the  computation  space,  was  successfully  used  to  model  the  CHAWS 
experiment.^  In  addition,  PIC  techniques  can  be  useful  when  developing  analytic  models.  In  order 
to  facilitate  these  modeling  techniques,  the  ability  to  perform  various  types  of  PIC  calculations  was 
built  into  Nascap-2k.  Recently,  additional  numeric  techniques  have  been  added  to  make  two  types 
of  PIC  calculations  more  flexible,  robust,  and  faster — Hybrid  PIC  and  Full  PIC.  In  a  Hybrid  PIC 
calculation,  the  problem  is  initialized  by  creating  ion  macroparticles  throughout  the  grid  to 
represent  a  constant  particle  density.  The  ion  macroparticles  are  tracked  for  one  timestep  and  then 
volume  potentials  are  computed  using  the  resulting  ion  density  and  a  barometric  electron  density. 

-((t)  +  e)/x^  for(t)>0 

-(eA^) 

exp((t)/6)  for(t)<0 
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where  <t)  and  0  are  the  plasma  potential  and  temperature  and  X  is  the  debye  length.  The  process  is 
repeated  for  the  time  period  of  interest.  In  a  Full  PIC  calculation  both  electron  and  ion 
macroparticles  are  tracked  and  volume  potentials  are  computed  using  the  resulting  plasma  density. 

Two  recent  enhancements  to  Nascap-lk's  PIC  capabilities  are  the  injection  of  macroparticles  from 
the  boundary  during  a  calculation  and  the  splitting  of  the  macroparticles. 

In  order  to  replace  macroparticles  that  are  collected  by  the  probe  or  escape  from  the  grid,  it  is 
necessary  to  periodically  inject  macroparticles  from  the  boundary.  This  allows  for  the  calculation  of 
current  for  longer  time  periods.  In  hybrid  PIC  calculations  without  boundary  injection,  the  low  field 
region  near  the  boundary  of  the  problem  develops  a  significant  negative  potential  due  to  ion 
depletion.  Boundary  injection  keeps  these  potentials  near  zero  by  replenishing  the  ions  that  have 
been  collected  or  escaped. 


Particle  injection  is  implemented  with  an  injection  point  at  each  quarter-boundary-surface-element. 
The  injected  particle  has  charge  equal  to  the  plasma  thermal  current  times  the  area  times  the  time 


interval  and  velocity  equal  to 


so  that  it  represents  the  inbound  half  of  the  plasma  density. 


When  the  spacecraft  is  moving  through  the  plasma  this  algorithm  is  modified  to  account  for  the 
motion.  The  current  and  velocity  are  computed  in  such  a  way  that  the  density  contribution  of  the 
injected  particles  varies  from  half  the  ion  density  for  a  stationary  object  to  the  full  ion  density  for  a 
high  Mach  number  object. 


Closely  connected  with  boundary  injection  is  macroparticle  splitting.  There  are  two  major  reasons 
for  splitting  macroparticles:  one  physical  and  one  numeric.  Even  at  moderate  potentials,  thermal 
effects  can  reduce  collected  currents.  Some  particles  near  the  sheath  edge  have  enough  thermal 
velocity  perpendicular  to  the  electric  field  that  angular  momentum  conservation  prohibits 
collection.  Particle  splitting  allows  for  a  representation  of  the  thermal  distribution  in  the  initial 
particle  distribution  and  in  particles  injected  from  the  boundary.  From  a  numeric  point  of  view, 
particle  splitting  can  be  used  to  keep  the  particle  weight  appropriate  to  the  grid  size  and  to  help 
maintain  the  smoothness  of  the  distribution.  A  large  particle  that  originated  in  an  outer  grid  is  split 
in  velocity  in  such  a  way  that  it  preserves  the  plasma  temperature  and  eventually  becomes 
distributed  over  several  volume  elements  in  an  inner  grid. 


Macroparticles  may  be  split  when  they  enter  a  more  finely  resolved  region  or  when  they  are  created 
either  at  the  boundary  or  throughout  the  volume  at  problem  initialization.  The  algorithm  used  is  as 
follows: 


1.  Particles  are  split  in  velocity  space  only.  Because  high-field  regions  are  often  of  interest,  spatial 
splitting  would  raise  problems  with  energy  conservation. 

2.  Each  particle  carries  a  temperature,  which  is  treated  as  isotropic.  The  fission  products  carry  half 
the  temperature  of  the  original  particle,  while  the  remaining  thermal  energy  appears  as  kinetic 
energy  of  the  split  particles.  That  each  macroparticle  has  a  temperature  means  that  they  can  be 
split  repeatedly,  each  time  the  particle  enters  a  more  finely  resolved  region. 

3.  For  splitting  purposes  the  Z-axis  is  defined  to  be  along  the  direction  of  the  particle  velocity,  the 
X-axis  randomly  chosen  in  the  plane  normal  to  Z,  and  the  Y-axis  mutually  perpendicular. 

4.  Particles  are  split  into  two  or  three  particles  along  each  axis,  except  that  a  particle  is  not  split 
along  the  Z-direction  if  the  kinetic  energy  exceeds  the  thermal  energy.  Not  splitting  along  Z 
helps  avoid  particle  proliferation,  but  makes  a  small  error  by  not  preserving  the  original  particle 
temperature  along  Z.  Eight,  nine,  or  twenty-seven  new  particles  result. 

5.  Particle  velocity  is  assumed  to  be  acquired  by  acceleration  rather  than  actual  drift  (i.e.,  spacecraft 
velocity).  If  there  is  actual  drift,  then  the  drift  velocity  is  removed  before  splitting  the  particle 
and  added  back  after. 
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6.  If  the  particle  is  split  by  two  along  the  X  or  Y  axis,  the  new  velocities  are  ±0.707^T/m  .  Along  the 

/  N 


Z  axis  the  velocity  increment  is  calculated  as  if  the  temperature  were  T-‘2muQ 


7.  If  the  particle  is  split  by  three  along  the  X  or  Y  axis,  there  is  a  zero-velocity  central  particle  and 
two  “probe”  particles  with  velocity  ±0.866-/f/m  .  Along  the  Z  axis  the  velocity  increment  is 


f 


\ 


calculated  as  if  the  temperature  wereT-2muQ 


The  new  particles  have  the  same  properties  as  the  original  particle  except  for  velocity,  weight 
(charge),  and  temperature. 


3  -  SELF-CONSISTENT  ANTENNA  CALCULATIONS 

We  performed  a  series  of  calculations  of  three-dimensional,  time-dependent,  self-consistent 
potentials  and  currents  for  a  spacecraft  with  a  VLF  transmitting  antenna.  The  Nascap-2k  model 
used  is  shown  in  Figure  1  and  Figure  2.  The  body  is  a  1.6  m  aluminum  cube  and  the  two  six-sided 
antennas  are  0.05  m  in  diameter  and  25  m  long  (3.75  m^  each).  The  grid  used  for  the  calculations  is 
shown  in  Figure  3  and  Figure  4.  The  mesh  unit  of  the  outermost  grid  is  2  m. 


Figure  1.  Nascap-2k  model  of  VLF  transmitting  spacecraft. 
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Figure  2.  Expanded  view  of  center  portion  of  Nascap-2k  model  of  VLF  transmitting  spacecraft. 


Figure  3.  Grid  used  for  VLF  antenna  calculations. 


Figure  4.  Close-up  of  center  of  grid  used  for  antenna  calculations. 


21 


5 


Results  from  four  calculations  are  shown  below.  The  parameters  of  the  calculations  are  shown  in 
Table  1.  In  all  cases,  the  potentials  are  adjusted  to  account  for  the  incident  current  each  timestep. 
The  incident  current  is  the  tracked  ion  current  together  with  an  analytic  electron  current  given  by 


1  = 


AJth  exp((})/0) 


Aith 


1  + 


4>' 


if  (j)  <  0 
if  (j)  >  0 


where  A  is  the  surface  area  and  jth  is  the  plasma  thermal  current.  One  antenna  is  floating  and  the 
other  has  a  variable  bias  with  respect  to  the  first.  The  aluminum  box  is  floating.  The  timestep  is  set 
so  that  macroparticles  move  a  reasonable  distance  each  timestep.  The  Hybrid  PIC  charge  density 
model  is  used. 

The  bias  applied  between  the  antenna  elements  is  the  sum  of  four  Fourier  components  that 
approximates  a  square  wave  with  amplitude  of  1  kV  and  the  indicated  frequency,  as  shown  in 
Figure  5.  In  the  absence  of  plasma,  each  antenna  element  switches  between  ±500  V.  The  square 
wave  excitation  is  desirable  in  order  to  minimize  time  spent  at  low  potentials,  for  which  the  plasma 
capacitance  is  highly  variable.  Having  a  phase-independent  capacitance  makes  it  easier  to  tune  the 
power  supply  for  optimal  operation. 


Table  1,  Parameters  used  in  antenna  calculations. 


C  use  1 

Case  2 

C  ase  3 

C  ase  4 

Density  (m'^) 

10* 

10* 

lo’ 

10’ 

Temperature  (eV) 

1 

1 

1 

IT" 

Species 

H" 

If 

IT 

If 

Ion  plasma  frequency  (kHz) 

2 

2 

6.6 

6.6 

Frequency  (kHz) 

10 

10 

12 

2 

Splitting  of  initial  macroparticles 

None 

All  grids 

All  grids 

All  grids 

Splitting  on  subgrid  entry 

None 

Three  levels 

Three  levels 

Three  levels 

Macroparticle  injection 

None 

Every  14 
timesteps 

Every  10  timesteps 

Every  timestep 

Number  of  macroparticles  (millions) 

0.50 

3.9  to  13.1 

3.9  tolO.4 

2.9  to  12.2 

Current  (mA) 

0.2 

0.03 

0.3 

0.25 

Phase  shift 

~0° 

10° 

12° 

-7° 
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Figure  5.  Applied  bias  values  and  resulting  antenna  potentials  in  the  absence  of  a  plasma. 


The  first  two  calculations  are  for  tenuous  plasma  (10^  m'^)  with  excitation  frequency  well  above  the 
ion  plasma  frequency.  Figure  6  shows  the  potential  of  the  antenna  elements  and  spacecraft  body  for 
Case  1,  which  has  no  particle  splitting.  In  the  initial  few  cycles  the  positive  antenna  element  collects 
substantial  electron  current,  driving  the  mean  antenna  potential  to  about  -500  V.  The  spacecraft 
body  initially  follows  this  transient  due  to  capacitive  coupling,  reaching  negative  potential  of  nearly 
200  V,  from  which  it  gradually  recovers.  During  the  switching  time  from  one  polarity  to  the  other 
there  is  a  period  during  which  both  antenna  elements  are  negative,  and  thus  collecting  some  ions  but 
no  electrons.  This  causes  the  positive  element  to  spike  at  positive  potential  following  the  switch, 
and  the  positive  potential  gradually  decays  during  the  half  cycle. 

Figure  7  shows  the  average  surface  current  to  an  antenna  element  for  the  case  with  no  particle 
splitting.  As  expected,  the  current  is  very  noisy  due  to  lack  of  splitting.  It  is  also,  at  0.2  mA 
(50  pA  m“^),  much  larger  than  the  2  pA  m"^  orbit  limited  current  for  a  cylinder  accounting  for  the 
thermal  angular  momentum.  This  is  as  expected,  since  without  splitting  the  thermal  velocity 
distribution  is  not  represented. 
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Figure  6.  Time  dependence  of  potentials  for  a  10*  m‘*  plasma  at  10  kHz  with  no  particle  splitting  (Case  1). 
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Figure  7.  Potential  and  collected  ion  current  of  antenna  element  2  for  a  10*  m’^  plasma  at  10  kHz  with  no 
particle  splitting  (Case  1). 

Case  2  has  the  same  physical  parameters  as  Case  1 ,  but  includes  particle  splitting,  thus  allowing  both 
for  smaller,  more  numerous  particles  and  a  representation  of  the  thermal  distribution.  Figure  8  shows 
the  potential  distribution  about  the  negative  antenna  element  when  fully  biased  at  2.1 1  ms,  at  the  end 
of  the  simulation.  At  this  low  plasma  density,  the  sheath  nearly  fills  the  grid  and  is  more  nearly 
spherical  than  cylindrical.  The  sheath  also  envelopes  the  positive  antenna,  providing  a  barrier  to 
electrons,  and  thus  allowing  the  system  to  float  more  positive  than  might  be  expected.  (This  barrier 
effect  on  electrons  is  not  modeled  in  the  present  simulation.)  The  collected  ion  current,  shown  in 
Figure  9,  is  smoother,  and  by  the  end  of  the  simulation  has  dropped  to  about  30  pA  (8  pA  m‘^),  which 
is  about  four  times  orbit  limited  for  a  cylinder.  This  is  reasonable  because  (1)  the  voltage  is  applied 
only  half  the  time,  and  (2)  additional  current  enters  through  the  “end”  of  the  sheath. 
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Potentials 


Figure  8.  Potential  structure  in  plane  through  center  of  spacecraft  and  antenna  elements  at  2.1 1  ms  for  a  10^ 
m*^  plasma  at  10  kHz  (Case  2). 

In  both  of  the  first  two  calculations,  the  plasma  current  to  the  negative  antenna  element  responds 
almost  immediately  to  the  switch  to  negative  potential,  even  though  the  applied  frequency  is  well 
above  the  ion  plasma  frequency.  This  occurs  because  the  applied  electric  fields  are  very  high,  and 
rapidly  accelerate  ions  to  speeds  far  in  excess  of  their  thermal  velocity.  However,  the  response  to 
the  positive  switch  is  less  sudden,  and  ions  continue  to  be  collected  by  the  positive  antenna  element 
through  nearly  all  its  positive  half-cycle.  The  fundamental  frequency  component  of  the  current  lags 
the  applied  voltage  by  a  10°  phase  shift,  substantially  less  than  the  ninety  degree  shift  expected  in 
the  linear  regime. 
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Figure  9.  Potential  and  collected  ion  current  of  antenna  element  2  for  a  10®  m'^  plasma  at  10  kHz  (Case  2). 

Cases  3  and  4,  at  high  and  low  frequency  with  respect  to  the  ion  plasma  frequency,  treat  a  plasma 
that  is  an  order  of  magnitude  more  dense  (10^  m"^)  than  the  first  two  calculations,  and  thus  has  a 
considerably  smaller  sheath  (Figure  10).  The  sheath  is  much  more  cylindrical  in  character  than  at 
lower  density  (Figure  8),  but  the  spherelike  end  cap  sheath  remains  substantial.  (Noise  in  the 
background  plasma  shown  in  Figure  10  is  at  the  IV  level.)  The  antenna  and  spacecraft  potentials 
(Figure  1 1)  follow  very  much  the  same  pattern  as  the  more  tenuous  plasma  (Figure  6),  except  that 
the  positive  potentials  are  better  held  in  check  by  the  higher  electron  current.  Similar  to  the  previous 
cases,  a  positive  spike  in  potential  is  seen  on  switching.  The  current  is  0.2  to  0.3  mA  (80  pA  m‘^), 
which  is  four  times  20  pA  m‘^,  the  orbit  limited  value  for  a  cylinder  in  this  plasma,  and  far  less  than 
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the  spherical  orbit  limited  value  of  600  pA  m’  .  The  current  variation  over  a  cycle  reaches  repeating 
on  the  fifth  cycle  at  12  kHz  and  on  the  third  cycle  at  2  kHz.  The  current  variation  with  time  (Figure 
12  and  Figure  13)  differs  between  the  high  and  low  frequency  cases.  Both  show  an  initial  rapid 
increase  as  ions  within  the  sheath  are  collected.  In  the  high  frequency  case  (Figure  12)  a  further 
gradual  increase  occurs  as  ions  near  the  sheath  edge,  which  take  longer  to  reach  the  antenna,  are 
collected.  The  fundamental  frequency  component  of  the  current  lags  the  applied  voltage  by  a  12° 
phase  shift,  more  than  that  seen  in  the  lower  density  case.  In  the  low  frequency  case  (Figure  13)  all 
ions  within  the  sheath  are  exhausted  early  in  the  pulse,  and  the  collected  current  is  limited  to  those 
ions  eroded  from  the  sheath  edge  during  the  pulse,  resulting  in  decreasing  current  during  most  of 
the  cycle.  Consequently,  the  fundamental  frequency  component  of  the  current  leads  the  applied 
voltage  by  a  phase  shift  of  about  7°. 


PotMTtialt 

Volts 
10, - 


30 


Figure  10.  Potential  structure  in  plane  through  center  of  spacecraft  and  antenna  arms  at  1.34  ms  for  the  high 
density  lO’ m’^  plasma  at  12  kHz  (Case  3). 
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Figure  11.  Time  dependence  of  potentials  for  the  high  density  10^  m‘^  plasma  at  12  kHz  (Case  3). 
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Figure  12.  Potential  and  collected  ion  current  of  antenna  element  2  for  the  high  density  10^  plasma  at  12 
kHz  (Case  3).  Note  that  the  current  lags  the  potential  by  a  12°  phase  shift. 
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Figure  13.  Potential  and  collected  ion  current  of  antenna  element  2  for  the  high  density  10^  m’^  plasma  at  2 
kHz  (Case  4).  Note  that  the  current  leads  the  potential  by  a  phase  shift  of  about  7®. 

Figure  14  illustrates  the  difference  between  the  high  and  low  frequency  behavior.  It  shows  the 
results  of  a  one-dimensional  calculation  of  a  spherical  antenna  in  a  plasma  of  density  4x10^  m”^ 
with  applied  sinusoidal  potential  of  5000V  at  10  kHz  (left)  and  1  kHz  (right).  At  high  frequency  the 
sheath  (defined  as  the  region  from  which  electrons  are  excluded)  extends  far  into  the  ambient 
plasma,  and  ions  are  collected  only  from  the  inner  region  of  the  sheath.  The  collected  current  can 
increase  during  the  time  the  potential  is  applied  as  ions  are  collected  from  progressively  farther  out 
in  the  sheath  at  each  cycle.  By  contrast,  at  low  frequency  the  sheath  region  is  nearly  emptied  at  each 
cycle,  so  that  the  current  must  fall  off  as  ions  from  the  sheath  edge  region  are  difficult  to  collect. 
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Figure  14.  One  dimensional  spherical  calculation  of  sheath  radius  and  ion  collection  radius  at  high 
frequency  (left)  and  low  frequency  (right). 


4  -  CONCLUSIONS 

We  have  used  Nascap-2k  to  study  the  plasma  interactions  of  a  high  voltage  near-square-wave  VLF 
antenna  in  MEO  plasma.  The  model  consists  of  two  25  meter  antenna  elements  biased  ±1000  V  at 
voltages  above  or  below  the  ion  plasma  frequency.  For  these  very  high  applied  voltages,  the  plasma 
response  is  nonlinear,  and  the  sheaths  are  more  spherical  than  cylindrical,  especially  at  the  lower 
density.  In  modeling  the  system  it  is  important  to  use  particle  splitting  and  injection  techniques  that 
replenish  depleted  plasma  at  boundaries,  maintain  appropriately  sized  macroparticles,  and  provide  a 
reasonable  representation  of  the  plasma  thermal  distribution.  When  this  is  done,  the  incident  current 
is  in  reasonable  agreement  with  orbit-limited  predictions. 

When  the  excitation  frequency  is  above  the  ion  plasma  frequency  the  plasma  current  to  an  antenna 
element  increases  sharply  when  the  potential  is  applied,  and  falls  off  relatively  slowly  when  it  is 
removed,  with  incident  ion  current  continuing  substantially  into  the  unbiased  half-cycle.  Below  the 
ion  plasma  frequency  the  ion  current  peaks  sharply  when  the  potential  is  applied,  falls  off 
substantially  during  the  half-cycle,  and  drops  sharply  to  zero  when  the  potential  is  removed.  The 
difference  can  be  understood  in  terms  of  the  sheath  being  depleted  of  ions  within  every  cycle  at  low 
frequency,  but  not  at  high  frequency.  The  current  lags  the  applied  voltage  for  excitation  frequencies 
above  the  ion  plasma  frequency  by  10°  at  the  low  density  and  12°  at  the  high  density,  which  is 
much  less  than  the  ninety  degree  shift  expected  in  the  linear  regime.  At  low  frequency  (calculated 
for  the  high  density  only)  the  current  leads  the  applied  voltage  by  7° 
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